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ABSTRACT 


The stability of the nonlinear dynamical system of two 
GRAVSAT-type satellites was investigated. The investigation was con- 
ducted by performing several numerical experiments which provide the 
simulations of the relative motions characteristics between the two 
satellites for various specified time intervals. The simulations 
include the relative range, range-rate, and the relative acceleration 
magnitude. 

These simulations were generated with respect to appropriate 
initial orbital elements which were obtained such that the instantaneous 
separation distance between the two satellites has small fluctuations 
from a specified constant separation distance. The simulation results 
indicate that the behavior of the relative motions is very sensitive to 
the initial orbital elements of the satellites and that for a specified 
time interval of interest, a stable behavior can only be achieved with 
the use of an appropriate set of initial orbital elements compatible 
with the gravity field used to derive them. 
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CHAPTER 1 


INTRODUCTION 

Over the past few yeari», observation of the relative motion 
between two satellites has gained a potential value for the study of our 
own planet. There are many geological and geophysical applications such 
as the study of the structure of the crust and the upper mantle, and 
accurate determination of the geold heights in the analysis of the 
Instantaneous and average shape of the sea surface [US National Academy 
of Sciences, 1979]. 

In the concept that is currently proposed, the accuracy of the 
coefficients which describe the gravity field of the earth could be 
improved by a simple analysis of the direct observations of the relative 
motions between two low altitude satellites separated by a specified 
distance in the same near-circular polar orbits. The relative motions 
which have been suggested for this analysis include range, speed and 
acceleration magnitude. In addition* the proposal would enable not only 
improved coefficients, but the determination of a more complete 
representation of the gravity field, including the anomalous potential. 

Tiie history of this concept is attributed to Wolff [ 1969 ], who 
proposed .a particular Satellite-to-Satellite Tracking (SST) configura- 
tion, known as the "low-low” system. In his proposal the relative speed 
along the line of sight between two satellites, as a first 
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approximation, la the difference in anomalous pot©i}tial sensed by the 
respective satellites. Therefore, as the Earth rotates, this pair would 
eventually cover the whole planet with observations of the gravity 
field. Since this original proposal, research efforts have been 
directed toward clarifying the theory behind SST, 

In the mid 1970’s, the idea of a dedicated gravitational 
satellite mission based on the SST principle began to develop and 
became known as the Gravity Satellite (GRAVSAT). In order to achieve 
the mission objective, it is essential that the satellites be made 
drag-free so that only gravitational forces act upon them. For this 
reason, a drag-compensatiQr technique has been proposed in which two 
proof-masses, one inside of each spacecraft, will be unaffected by sur- 
face forces through periodic activation of thrusters which would offset 
the sensed surface forces acting on the spacecraft enclosing the proof- 
mass. Such a system has been constructed by Stanford University, 
referred to as the Disturbance Compensation System (DISCOS) [Pisacane et 
al., 1981]. This drag-compensation system of satellites would be capa- 
ble of remaining in low altitude orbit for a nominal mission lifetime of 
up to six months. As long as the drag-compensation mechanism is oper- 
able, the atmospheric resistance forces will be effectively eliminated 
and the satellites will move in orbits governed by gravitational forces. 
Without such a mechanism, the large drag forces would quickly decay the 
orbits and make the desired determination of the gravity field far more 
complicated. 
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Tfie line of sight measurements between the two satellites 
would be made by an extremely accurate radar interferometric technique 
developed at the Applied Physics Laboratory of John Hopkins University 
[Pisacane et al., 1981 Ihe range-rate precision is expected to be on 
the order of ym/seo between the two satellites, 

The general low-low configuration system along with the 
desired mission characteristics are shown in Figure 1. The proposed 
configuration includes a separation of 300 km between the satellites at 
an altitude of 160 km# Although it is desirable in some cases for the 
respective satellite orbit planes to coincide, they are illustrated as 
noncoincident for the general case. 

In the following section, descriptions of the major approaches 
for some of the recent studies, conducted under the use of simplifying 
assunptions, are presented. These studies provide the motivation for 
this investigation. 

1. 1 Motivation for this Study 

Basically, an analysis of data for the improvement of geopo- 
tential knowledge consists of three major sections, which can be stated 
as follows: 

1. Develop observation equations in order to relate the measurements 
to the anomalous potential of the Earth. 




ii . Process the data . 


iii. Determine the accuracy of the results. 

Of course, the accuracy of the observations and the overall modeling 
will determine the accuracy of the results. 

In general , any globally consistent solution for a degree and 
order 180 gravity field in terms of spherical harmonic coefficients by 
least-squares techniques must solve for approximately 32,000 parameters. 
Consequently, features of the gravity field with 100-200 km in 
wavelength can be determined [Pisacane and Yionoulis, 19801, To avoid 
the massive matrix inversion computations, associated with such a large 
nunber of parameters, different methods have been proposed. For exam- 
ple, Colombo [1981] takes advantage of the orthogonality properties with 
respect to latitude of the spherical harmonics. Also, Kaula [1982] pro- 
posed the use of Fourier analysis for each revolution of the satellites. 

Some common simplifying assumptions have been made in studies 
by Colombo [1981], Gaposchkin [1980] and Kaula [1982]. A collection of 
these assumptions are summarized below; 

1. Both satellites describe coplanar, circular , polar orbits with the 
same geocentric radius. 

2. Ihe satellites are separated by a constant distance and the separa- 
tion distance has small fluctuations. 
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3. The plane of the orbits is fixed in inertial space, and fluctua- 
tions in the geocentric angle between the two satellites are small. 

These simplifying assumptions lead one to ignore the actual perturba- 
tions of the orbits in the radial, transverse, and normal directions 
which are produced due to the effects of the anomalous potential. 

The study of the validity of these assumptions and to deter- 
mine the extent and rate at which the satellites might drift apart over 
various specified time intervals was the primary objective of this 
investigation . 

1,2 Description of the Present Study 

In this study, attention has been focused on the relative 
motion of the satellites, and basically, the stability of this dynamic 
system as far as the separation distance was concerned. To support The 
objective of this study, a set of initial orbital elements for each of 
the satellites was found so that the instantaneous separation distance 
(or geocentric separation angle) has small fluctuations from a specified 
constant separation distance. Therefore, a clear prospective of whether 
the stated approximations are valid can be concluded. A brief descrip- 
tion of the material presented in this investigation is given as fol- 
lows . 

Chapter 2 presents the derivation of the equations of motion 
considering only the gravitational forces acting upon the satellites. 
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In addition, a rotating coordinate -s^ystem attached to the trailing 
satellite was chosen so that the relative motion between the two satel- 
lites can be discussed qualitatively. 

Chapter 3 provides the mathematical formulation of an optimi- 
zation technique used to solve for the initial orbital elements of the 
satellites. In this method the initial orbital elements of the trailing 
satellite are specified. Then the orbital elements of the lead satel- 
lite are determined such that the relative range between the two satel- 
lites is nearly constant over a specified time interval. 

Chapter 4 includes the numerical simulations of the relative 
motions between the two satellites. 

Finally, in Chapter 5, based upon the analysis of the nuneri- 
cal simulations, the conclusions of this investigation are summarized. 



CHAPTER 2 


DYNAMICAL MODEL 

This chapter presents the differential equations which 
describe the motion of two Satellites orbiting the Earth under the sole 
influence of the earth’s gravitational field. Section 2.1 includes the 
gravity field models used in this study. Section 2.2 provides the rela- 
tionships between the different coordinate systems used in this investi- 
gation and the equations of motion are presented in Section 2.3. Sec- 
tion 2.4 provides an analytical definition of the relative range-rate 
between the two satellites and a satellite-centered coordinate system is 
Introduced. In Section 2.5t an observation equation which relates the 
line of sight acceleration magnitude and the gravity potentials is 
obtained such that no simplifying assumptions have been used. 

2. 1 The Geopotential Models 

The gravitational field of the Earth as measured from above 
its surface can be represented by a potential, U, which satisfies 
Laplace's equation; 

U = 0 (2.1) 

The solution to Eq. (2.1), using the spherical harmonics 
representation, is given by the following infinite series [Kaula, 1966] 
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assuming an origin located at the center of mass: 


U = 


GM 


N 

1 - 2 

n=2 


Pn(3in(|.) 


( 2 . 2 ) 


N n 
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Universal gravitational constant 

mass of the Earth 

equatorial radius of the Earth 

radial distance from the center of the Earth 

geocentric latitude with respect to the equ^^tor 

angle of longitude with respect to the Greenwich 
meridian, measured eastward 

Legendre polynomials 

associated Legendre functions 

spherical harmonic coefficients 

order of the spherical harmonics 
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m s degree of the spherical harmonics 

The first term on the right hand side of Eq. (2,2) represents 
the central body potential and the remaining terms are due to the non» 
sphericity of the Earth. 

Eq. (2,2) is an approximate model of the geopotential where N 
denotes the truncation order of the geopotential model. For most 
approximations, N is some finite Integer. The spherical harmonic coef- 
ficients have been estimated from various data sets. A brief descrip- 
tion of some recent determinations are as follows: 

• Goddard Earth Model (GEM) 9 CLerch, et al . , 19793 is a gravity 

model based solely on optical, laser, and electronic observations 
taken on 30 satellites. It has harmonics complete to degree and 
order 20 (20x20) with selected higher degree terms. 

• GEM-10 CLerch, et al.» 19793 is a combination solution containing 
a global set of 5 degree surface gravity anomalies along with the 
data in GEM-9. It has harmonics complete to 22x22 with selected 
higher degree terms. 

• GEM-10B CLerch, et al . , 19813 is a gravity model that combines the 
GEM-10 data set with those obtained from altimeter data of the 
GEOS-3 geodetic satellite. GEM-1 OB represents a significant 
improvement over GEM-9 and GEM-10. It consists of a field com- 
plete to 36 x 36 in harmonics. 
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Among the above gravity models, the GEM-9 and GEM-10B were used in the 
nunerical analysis of this investigation. 

2.2 Coordinate Systems 

Several different coordinate systems were used in the 
mathematical formulation of the dynamical system for this investigation. 
Figure 2.1 illustrates these coordinate systems. For simplicity, the 
precession, nutation, polar motion and changes in rotational speed of 
the earth have been ignored. The coordinate systems are defined as fol- 
lows . 

The nonrotating, inertially fixed rectangular coordinate sys- 
tem with the origin at the Earth' s center of mass denoted by {X,Y,Z} , 
has the X-axis directed along a fixed direction (e.g., the mean Equinox) 
and the Y-axis is in the equator. All of the differential equations 
which describe the motion of the satellites are expressed in this iner- 
tial reference system . 

A rotating, Earth-fixed rectangular coordinate system with the 
origin at the Earth's center of mass is present to account for the 
Earth's rotation and the description of the gravity field coefficients. 
This system, denoted by {x,y,z}, is a right-handed system with x-y plane 
coinciding with the equatorial plane of the earth. Therefore, it 
rotates about the Z-axis with the angular velocity of the Earth. The 
magnitude of this rotation angle at any instant can be computed by 
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a r .. 0)^^ ( t - ) 


(2.3) 


where is the magnitude of the Earth's rotation about, the Z-axis, and 
cx^ is the initial angle at the time t^. The angle a is measured east- 
ward from the X-axis to the Greenwich meridian. Thus, the components of 
any vector in the inertial coordinate system can be transformed into the 
Earth-fixed system by the following transformation matrix: 


T = 


cos a sina o 

-sina cosa 0 

0 0 1 


(2,iO 


Of course, T is an orthogonal matrix, therefore the relation 



can be used for the inverse transformation. 


Finally, the gravitational field of the Earth as defined in 
Eq. (2.2) is represented in an Earth-centered spherical coordinate sys- 
tem denoted by {r,({).X}. The coordinates of any point on a satellite's 
orbit in this reference frame is related to the Earth-fixed rotating 
coordinate system by the following relations: 


X s r cos4) cosX 

y = r coS(f) sinX (2.5) 

z = r sinc|) 

where r is the radial distance from the center of mass of the Earth, tj> 
is the angle of latitude from the equator, and X is the angle of longi- 
tude measured eastward from the Greenwich meridian . 
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In addition to the above coordinate systems, a rotating refer- 
ence frame attached to the trailing satellite is introduced in Section 
(2.4). This system is used for the theoretical discussion of the rela- 
tive motion between the two satellites. 

2.3 Equations of Motion 

In this section, differential equation which describe the 
relative motion of the two satellites considering only the gravitational 
forces are presented. 

Assume that the mutual gravitational attraction between the 
two satellites is negligible compared with the earth’s gravitational 
forces acting on them. Therefore, considering the satellites as point 
masses, the total inertial accelerations of the satellites subject to 
the Earth’s gravitational field are given by 
• • 

R^(t) = V U^(t) (2.6) 

• • 

^^(t) = V 

where the numerical subscripts 1 and 2 denote the trailing and the lead 

satellite respectively; ^ is the acceleration vector; ^ is the gradient 
operator; and U is the gravitational potential of the Earth defined by 
Eq. (2.2). 

Now, with the use of appropriate transformations, the carte- 

I 

Sian components of the acceleration of each of the satellites can be 




obtained by the following equations: 


(2.8a) 


(2.8b) 


(2.8c) 


The three second-order differential Equations (2.8a) through (2.8c) can 
be converted to six first-order differential equations as follows: 


l?(t) = V(t) 


(2.9) 


V(t) =7 U(t) 


( 2 . 10 ) 


where R(t) and V(t) are the position and velocity vectors, respectively. 

Given a set of initial components of position (X^.Y^.Z^) and velocity 

• • • 

(X_,Y_,Z_) vectors for each of the satellites, Equations (2.9) and 

Q Q Q 

(2.10) can be nunerically integrated to provide the rectangular com- 

• » e 

ponents of position (X,Y,Z), and velocity (X,Y,Z) vectors at any time 
for each respective satellite. Ihen, the relative range vector (T) can 
be obtained by 


p = ^2 ~ 


( 2 . 11 ) 


where the subscripts 1 and 2 denote the trailing and the lead satellite 
respectively. So , the relative range between the two satellites (p) can 
be computed by 




( 2 . 12 ) 
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Dlffereiitlatlng Eq. (2.11) twloe with respect to time and 
using Eqs. (2.6) and (2.7) yields 

• • 

p = 7 Ug - V (2.13) 

• • 

whore p is the relative acceleration between the two satellites, 
expressed in the inertial system. Eq. (2.13) is a nonlinear differen- 
tial equation which describes the relative motion between the two satel- 
lites. It should be noted that if the initial orbital elements of the 
satellites are not chosen properly, an irregularly secular behavior in 
the separation distance between the two satellites could result, 

2.H Definition of the Relative Range-Rate 

The relative velocity vector (p is the time derivative of the 
relative range vector. Iherefore, differentiating Eq. (2.11) yields 

• « • 

p = T?2 - (2.14) 

in the inertial system. 

• ' 

The relative range-rate between the two satellites (p) is 
mathematically defined as the projection of the relative velocity along 
the line of sight between the two satellites. In other words, it is the 
time rate of change in the distance between the two satellites. Accord- 
ingly, 

P = P • P = ^2*15) 

IIpII 
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where p is a unit vector along the line of sight, and ))p|| is the rela- 
tive range given by Eq, (2.12). Tlius, it is of importance to note that 

p = |jpj| but p is not the same as the magnitude of p . Now, by substi- 
tution of Eqs. (2.12) and (2.14) in Eq. (2.15)» and using the rectangu- 

lar components of p in the inertial coordinate system* the range-rate 
can be computed from 

p : (2,16) 

[(Xj-x,)® * 

Assume a relative reference frame is attached to the trailing 
satellite with one axis directed along the line of sight toward the 
leading satellite as shown in Figure 2.2 (more detail on this reference 
frame will be given in Section 2.5)* then the relative velocity vector 

9 

in this reference system (p ) can be obtained by 

r 

p^ = p - si X p (2.17) 

r 

where l2 is the total angular velocity of the satellite-centered rotating 

coordinate system, and p is the relative velocity vector expressed in 
the inertial reference frame given by Eq. (2.14). Taking the dot pro- 
duct of p viith Eq. (2.17) results in the following equality: 


P * Pr " P*P “ P 


(2.18) 


since 


original Pm : m 

OF POOR QUAUfy 


Figure 2.2 


Satellite-Centered Orthogonal Coordinate System 


Trailing Satellite 


Lead Satellite 
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p • X p) = 5 • (p X p) = 0 . 

Nots that Eq. (2.18) yields again the relative range-rate. Therefore, 
it can be concluded that the relative range-rate can be defined as the 
projection of the relative velocity in either the inertial or the rotat- 
ing reference frame along the line of sight. 

2.5 Definition of the Relative Acceleration Magnitude 

In the following discussion an observation equation is 
presented in which relative range, relative range-rate, and the line of 
sight relative acceleration are directly related to the inertial com- 
ponents of the position, velocity, and the acceleration vectors of the 
satellites. This equation was originally developed by Shao [1982]. In 
deriving this equation, it should be noted that the line of sight meas- 
ured quantities are different from the inertial components of the vec- 
tors. This equation is of significance due to its generality and that 
no simplifying assumptions or approximations have been made. In addi- 
tion, the coordinate system used in the derivation of this equation is 
identical to the one described in Section 2.4, This choice of reference 
frame is important because it is inherent to the dynamics of the system, 
and furthermore the rotation of the lead satellite about the trailing 
satellite is described by the angular velocity of this reference frame. 
Further details about the properties of this relative reference frame 
are given in the following discussion. 
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■Rie unit vector p is defined to be a unit vector along the 
line of sight, so its time rate of change can be obtained from 

• • 

p = I? X p 4. 3^ . (2.19) 

Since p is a unit vector, its time derivative in the rotating frame does 
not change; in other words, p^ = 0. Thus Eq. (2.19) reduces to 

e 

p = X p . (2.20) 

However, a constant length vector is perpendicular to its time deriva- 
tive. Therefore, 

» 

3.3=0 (2.21) 

In addition, by the use of vector algebra equalities, it can be shown 
that p is perpendicular to fi. The proof is as follows: 

p*J2=CfiXp).j^ 

= p • ^ ^ Q ) 

Since x =0» the above equation reduces to 

3 . si = 0. (2.22) 

Equations (2.20)* (2.21), and (2.22) represent necessary and sufficient 

conditions for 3 * p» to be mutually perpendicular; conse- 

quently, they form an orthogonal basis for the previously described 
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t » • .j 




as.' 


satellite-centered reference frame. 


Recall that the relative motion between the two satellites was 
given by Eq. (2. 13) or 

p = V Ug - V = h-h 

where and are the potential functions evaluated at the respective 

satellite positions. Note that the left hand side of Eq , (2.23) is the 
relative acceleration vector in the inertial coordinate system, and it 
can also be expressed in the relative reference frame as follows; 

• • • 

d , -r ^ 

P 5 ( P ) 

= -JE ( Pr ♦ fj X P ) 

• • • • • 

= (p^+fixp^) + (fixp + sTxp) 



Substitute + Ti X p for p in the term p x p, then it follows that 

p = pp + 2( X p^ ) - p + n X p (2.2U) 

Forming the dot product with ^ and both sides of Eq. (2.24) yields 

p . p = p . p^ + 2 p . ( D X "p^ ) (2.25) 

”P*(f2 p) + p*(f2xp) 

Note that each term on the right hand side of Eq. (2.25) simplifies to 


pn 
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P . P. 


2 p . ( Q X ) s 2 Q . ( X p ) = 


A / 2 \ 

p • ^ p ^ 


2 

n p 




= 0 


(2.26a) 

(2.26b) 

(2.26c) 

(2.26d) 


Substitution of Equations (2.26a) through (2.26d) in Eq. (2.25) results 
in 



(2.27) 


Using Eq. (2.23) and also rearranging the terms, Eq. (2.27) can be writ 
ten in the form 


p = p + ( Rg " ^1 ^ • P 

or 

p = p + ( V Ug - V ) . p (2.28) 

Eq. (2.28) is the observation equation v^ere p is the line of sight 
relative acceleration, p is the relative range, and is the magnitude 
of the angular velocity of the rotating relative reference frame, which 
can be computed in terms of inertial components of position and velocity 


of the two satellites as follows. 
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Taking the cross product of p with both sides of Eq. (2.20) 

yields 

p^p”p^Cj2*p)* (2.29) 

Using the following vector algebra equality 

Xx(lx'S) = F(7[.‘(J)~7(7[.¥) 
then Eq. (2.29) oan be written as 

pXp=Q(3x3)-p(3.fi). (2.30) 

Since p is perpendicular to Eq. (2.30) reduces to 

p X 3 = . (2.31) 

— 0 
Therefore, is defined by Eq . (2.31). and fj- can be computed from 

= II P « p 11^ 

• • 

= ( p X p ) . ( 3 X p ) (2.32) 

Using the following vector algebra equality 

( Ax'S ) .( 1CxT5 ) s ( Hl.T! )( B.'D ) - ( l.D )( 'H.'C ) 

then Eq. (2,32) can be rewritten in the form of 

= < p.p H p.p ) - ( p.p )( p.p ) (2.33) 


< 0 - 
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Note that p is perpendicular to p, thus p.p = 0. Then, Eq. (2.33) 
reduces to 


• • 


—2 _ A * _ / a \2 

SI = p . p s Cp; 


(2.34) 


is the time rate of change of p, and can be obtained as follows 


^ _ U / A \ 

p " ■ar ^ p ^ 


d (~ ) 


p p - p p 


(2.35) 


Finally, can be computed by substitution of j5 in Eq. (2.34), as fol- 
lows ; 




12 


P P - P P 




,4 L 


p2 ( p )2 _ 2 P^ P^ + P^ P^ 


(Q)^ - P^ j 


(2.36) 


where P is the relative range-rate given by Eq. (2.16), and 


_p • *P ’'o • *? 

(p)^ = (Xg - x^)"^ + (Yg - Y^)'^ + (Zg ‘ z^r . 
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Note that in general is time dependent and a function of 
position and velocity of the satellites. And in particular, if the 
coplanarity, circularity, and constant relative separation orbital 
assumptions are made, then Tl would be constant with a fixed direction 
perpendicular to the orbital plane. 


Now, using Eq. (2.36) for substitution of in Eq. (2.28), 

• • 

then P can be obtained in terms of the rectangular components of the 
position, velocity, and acceleration of the satellites in the fixed 
coordinate system as follows: 


t [ 


(P)^ - P^ 


‘2 

p) - P 


+ ( V U2 - V ) . p 


. i[ 

’ (P)^ - p^ + p . p j 

= i [ ^^2-^ 


+ ( p ) . p 


X^)2 + (Yg-Y^)^ + ^ 


(2.37a) 


(2.37b) 


(2.37c) 


(2.37d) 


• • • • 


• • • • 


+ (X2-X^)(X2"X^) + (Yg-Y^XYg-Y^) + (Z2' 


-z,k‘4:z,) 1 


From any of the above equations, it can be concluded that the 
line of sight relative acceleration is different from the magnitude of 
the relative inertial acceleration. Furthermore, their difference is 
not a constant quantity. 
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In summary, the set of equations required to calculate the 
relative motions between the two satellites are as follows. The state 
of each of the satellites is calculated by integration of Eqs. (2.9) and 
(2.10) for the respective satellite. Then, the relative range is com- 
puted by Eq. (2.12), the relative range-rate by Eq. (2.16), and the 
relative acceleration magnitude by Eq, (2.37d). 



CHAPTER 3 


DETERMINATION OF APPROPRIATE INITIAL ORBITAL ELEMENTS 

In the previous chapter, it was shown that the relative motion 
between the two satellites is governed by the nonlinear accelerations 
produced due to the earth’s gravitational field. 

The stability of the present nonlinear dynamical system has a 
direct dependency on the initial orbital elements of the satellites. As 
it will shown in Chapter 4, a simple nunerical simulation of the rela- 
tive range between the two satellites, with any arbitrarily selected 
initial conditions, will illustrate a secular behavior in the separation 
distance with oscillations of large amplitude. Therefore, in order to 
obtain a stable behavior in the relative motion, and to maintain small 
fluctuations between the two satellites, a method is developed to obtain 
appropriate initial orbital elements. In this method, the initial orbi- 
tal elements of one of the satellites is specified. Then, using the 
principle of least squares, the initial orbital elements of the other 
satellite are nunerically obtained to minimize the sum of the squares of 
the changes in the separation distance from a desired constant distance. 
The formulation of this problem is presented in Section 3.I, and the 
solution is obtained by using an iterative method which is introduced in 
Section 3*2. Finally, a computational algorithm is summarized in Sec- 
tion 3 . 3 . 
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3. 1 Problem Formulation 

Having described the dynamical system represented by the sys- 
tem of Eqs. (2.6) and (2.7), the differential equations describing the 


state (position and velocity) 

of the satellites, X^ and X^, can be writ- 

ten in the 

form of 


X^(t) 

= F, (X, ,t) 

(3.1) 

• 

. Fj<Xj.t) 

(3.2) 


where F^(.) and nonlinear functions of and respectively. 

Let the initial components of position and velocity of one of 

the satellites, say the trailing satellite, be fixed and denoted by X^ . 

o 

Then, it is desired to solve for the appropriate initial components of 
the position and velocity of the lead satellite as arranged in the esti- 
mation state vector 



Let be the desired constant separation distance between the 
two satellites. 'Ihen, at any Instant of time t during a specified time 
interval, the change in the separation distance is defined by 

r = P - (3.4) 
where p is the relative range which can be computed by Eq. (2.12) and r 
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is referred to as the residual . 


The optimal solution is obtained by solving a deterministic 
optimization problem for which the performance index is the sum of the 
square of the residuals. In this approach, the commonly defined perfor- 
mance index can be expressed as 


P = t r^ } dt 
t 

o 

where 


P = the function to be minimized 


(3.5) 


t = the initial time 
0 

= the specified final time 


First variations of P with respect to the components of the 
estimation state vector are obtained by 



"f 

S t 

t 



r } dt 


(3.6) 


In order to minimize P, first variations of P must vanish. Hence, six 
nonlinear equations are produced which are 



8p 



r } 

J 


dt 



3p 



r } dt 



(3.7a) 


(3.7b) 


i* 1 ■». 
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‘'f 

S ( 

t 


ap 

T2T 


r } dt = 0 


X ( 

t. 


8p 

"• i'"“ 

3Xo 


r } dt = 0 


X i 

t 

0 


3p 


r } dt = 0 


X t 

t. 


3p 

3 Z 


r ) dt = 0 


Equations ( 3 . 7a) through (3.7f) can be put in a matrix form as 


f = 


‘'f 

X t 

t. 


ap 


3X, 


r ) dt = 0 


and If the following definition is introduced, 


H = 


3p 


3 X, 


Eq. (3.8) then can be expressed in terms of H matrix by 


f = X t = 0 
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(3.7c) 

(3.7d) 

(3.7e) 

(3.7f) 
follows ; 

(3.8) 

¥ 

(3.9) 

( 3 . 10 ) 
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Thus, the problem reduces to solving the six nonlinear equations of T in 

terms of six unknown components of the estimation state vector of JCg • 

‘^0 


3.2 Solution Method 


A very frequent computational problem is to find some or all 
of the solutions of a system of n simultaneous nonlinear equations. One 
approach to solve such system of equations is to generalize to n dimen- 
sions one of the iterative processes used for solving a single equation; 
see for instance Forsythe, et al. [197?]. There are several methods 
available from which the system of equations defined by Eq. (3*10) can 
be solved. However, the ttewton-Raphson iterative method was chosen and 
the solution is obtained as follows. 


Assume that the partial derivatives of the functions f^^ with 

respect to the components of X, can all be computed. Then, let J(X« ) 

'^o 

represent the Jacobian matrix whose (i,j)th element is defined by 






(3.11) 


As in the case of one dimension, the Newton idea is to start with an 

arbitrary Y , say "X? . Then, the function f can be linearized at "X? 

^o ^o ^0 

by expanding 7 into a Taylor series which by using the above definition 
it can be written in the form of 


fCXg ) = f(x| ) + J(X| ) (Xg - X^ ) + . . . 

“^0 o 000 


(3.12) 
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Now, keeping only the terms of degree 0 and 1, the first approximation 
to the solution of Eq. (3.10) can be obtained by solving the following 
equation for : 


T(X| ) + J(X^ ) (Xp - ^ ) - ^ (3.13) 

and the first approximate solution s is obtained which is: 

o 

Xg = ^ ) f(X^ ) 

0 0 0 ^o 

Of course, in the general step of the iteration, if Xp is the solution 

■^o 

th 

obtained at the k iteration, the next approximation to the solution 
can be computed by 


;k+1 


? x: 




vk 

2. 


f(Xp ) 
0 


(3.14) 


Note that the second term on the right hand side of Eq. (3.14) 
represents the correction to the approximate initial state of the lead 

satellite. And if the correction at each iteration is denoted by <S(Xp ) 

then Eq. (3.14) can be written in the form of: 


7k+1 

^0 


Xg - S (Xg ) 
o o 


where 


6 (Xp ) 


j’^Xp ) 


f(Xp ) 
rk ^o 


(3.15) 


(3.16) 


X*^ 

2 . 
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Existence of the Solution 

Convergence of the solution using Newton's method depends on 
the initial guess of the state vector which is required to start the 

— • 

iteration process. In other wrds, suppose x is a nominal solution of 

^o 

f s 0. Then, if the initial guess ) is sufficiently "close" to ^ , 

*^0 ' p 

the Newton iteration will converge, and 
Xp -> }^ as k -> 00 

Moreover, high-speed convergence will start when gets sufficiently 

^0 

close to X* . Although, it should be noted that as soon as the differ- 

t h 

ence between the obtained solution at the k iteration and the nominal 
solution defined by 



reaches to the order of the distance between nearby floating point 
numbers, the granular structure of the number system makes it no longer 
possible to continue the iteration. 

The other criteria which should also be satisfied at each 
iteration is that the Jacobian matrix be non-singular so that it can be 

converted in order to compute the correction ). 
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Computation of the Jacobian Matrix (J) 
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Ihe Jacobian matrix (6x6) was defined by Eq. (3.11) and can be 
computed as follows. Using Eq. (3.8) for substitution of T into Eq. 
(3.11). it then follows: 


J(Xp ) 


3 X. 


tf 

S ( 

t^ 


3p 

9X, 


r }dt 


tf 

s {— 


3p 


3 In 


r } dt 


(3.17a) 


(3.17b) 


''f 

^ S t 

t. 


8p 


8p 

Big 


Big 

0 


0 


} dt 


Neglecting the small terras (second partial derivatives of P with respect 
to the components of X ) which are associated with the first terra on 

the right hand side of Eq. (3.17b), and also using the definition of Eq. 
(3.9). J(Xn ) can then be coraputed approximately by 


JOT ) Z J { H } dt 


(3.18) 
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J is a symmetric matrix and along with T in Eq, <3.10), both 
are expressed in terms of H which can be obtained as in the following. 


Computation of H Matrix 

TTie components of the H matrix (1x6) are not an explicit func- 
tion of the components of . Therefore, using the chain rule, the H 

o 

matrix can be written in the form of 


8p 


3p 


DXp 

SXj 


axg 


axj 

O 


- 


0 


And, if the following definitions are introduced. 


(3.19) 



( 3 . 20 ) 


$(t,t ) 
0 



aXg 

o 


( 3 . 21 ) 


then Eq. (3.19) can be expressed in terms of H, and $ as 

H = H $ (3.22) 

Components of the H matrix can be obtained by taking the par- 
tial derivatives of p, defined by Eq, (2.12), with respect to the com- 
ponents of the state of the lead satellite ^2 at some instance of time. 
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These partials are: 

8p ^2-X 
= — 

3P ^2 
3Y2 - p 

80 hjlL 

3Z2 " P 


1 


1 


1 


8X2 

,^Ji s 0 



(3.23a) 

(3.23b) 

(3.23c) 

(3.23d) 

(3.23c) 

(3.23e) 


Collecting the partials from Eqs . (3.23a) through (3.23e), the H matrix 
is then computed by 


H = 


V^1 

P P P 


0 


The State Transition Matrix 


(3.211) 


The state 
Based on Eq. (3.2) 

satellite, assuming 


transition matrix , is defined by Eq. (3.21). 
which represented the nonlinear dynamics of the lead 

that Xg(t) is a nominal solution of Eq. (3.2), then 


the 




deviation of from the ’’true” solution Xg(t) at any instance of 


time is defined by 
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= XgCt) - X*(t) 


Furthermore , assuming that X 2 (t) is sufficiently "close” to XgCt) , then 


— _» _ 
the function Fg can be linearized about Xg by expanding Fg into a Taylor 

series. Keeping only the terms of degree 0 and 1, and provided that 

A(t) matrix is obtained from Fg(.) by the following definition; 


A( t) 





(3.25) 


then it can be shown that the following linear system of equations are 
satisfied [Brogan, 197^3: 

7g(t) = A(t) XgCt) (3.26) 

and that the state transition matrix satisfies the following differen- 
tial equations: 

<f(t,t^) S A(t) $(t,t^) , $(t^,t^) S I (3.27) 

t^ere I is an identity matrix. Thus, the state transition matrix can be 
computed by numerical integration of Eq , (3.27) from t^ up to t. 

Summarizing, in order to minimize the performance index 
defined by Eq. (3.5), the following differential equations need to be 
simultaneously Integrated to result in the appropriate initial orbital 
elements of the lead satellite: 
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• Six equations of motion for the trailing satellite given by 


. Fl . X^(V = X^ 


(3.28) 


• Six equations of motion for the lead satellite given by 


h ‘ ^2 ' = *2 <3.29) 

o 

• Six equations for the elements of 6x1 function T given by 

f = r , f(t^) = 0 (3.30) 

• Thirty-six equations for the elements of the 6x6 Jacobian matrix 
given by 

j = H , J(t^) = 0 (3.31) 

• Thirty-six equations for the elements of the 6x6 state transition 
matrix given by 

= A <^ , $ (t^,b^) S I (3.32) 

•One optional equation may also be integrated to determine the 
effect of the convergence of the solution on the value of the per- 
formance index which can be computed by integrating the following 
differential equation; 

P = r^ , P(t^) =0 (3.33) 
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3 . 3 Computation AlRorithm 

In this section a batch type computational algorithm [Tapley, 
1973] is presented which requires the integration of 91 first-order dif- 
ferential equations summarized by Eqs. (3.28) through (3.33). The fol- 
lowing steps lead to the iteration process of the estimation of the 

appropriate initial orbital elements for the lead satellite, X* . 

0 

0» Initialize the state vectors, X, (t ) = X, and X«(t ) s xj , Set o 

lo 1 2 0 2 ^c 

o o 

to the desired constant separation distance. 

1. Set the iteration nunber ks1 and X^ s X° . 

o *^0 

2 . Initialize the matrices, (6x1) "fsO, (6x6) JsO, (6x6) $=1, and PsO. 

3. Set is1 and t^ = t^. 

4. Compute the A matrix at t^^ . 

5. Compute the H matrix at t^ . 

6. Compute the H = tf *5. 

7. Compute r = tl^(t.) - X^( tp j ) - p^. 

8. Increment time: t^ = t^ + At. 

9. Check: if t^^ > t^, then go to Step 10. Otherwise, 
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• Integrate the equations of motion to get t and • 

• Integrate the state transition matrix to get at t^ . 

• Integrate T = r . 


• Integrate J = H H. 


• Integrate the performance index, P = r*^ (optional). 


• Go to Step 4. 


Tk 

10. Compute the correction o at t^, 


5 = J"^ f. 


IT. Compute 1 ^*'^ = Xg 
o o 


12. Check the convergence. If > tolerance, then; 


• Increment k: k s k+1. 


**k ^k4* 1 

• Replace X- by Xp , go to Step 2 for the next iteration. 


13. Otherwise, the final approximate solution is obtained which is 

o 

the same as Kp"*"^ . 

‘^0 


CHAPTER ^ 


SIMULATIONS OF THE RELATIVE MOTIONS 

In this chapter f the relative motion characteristics between 
two GRAVSAT type satellites have been numerically investigated. Tlie 
investigation includes the analysis of the changes in the relative 
motions produced due to the effects of the earth’s gravity field models 

I 

which were described in Section 2.1. 

Several numerical experiments were performed in order to pro- 
vide the simulations of the relative motions which are the relative 
range, range-rate, and acceleration magnitude. The equations used to 
compute the relative motions were derived in Chapter 2, where they are 
given by Eqs. (2.12), (2.16), and (2.37d), respectively. 

The estimation algoritlvn which was developed in Section 3.3 
has been used extensively to obtain the appropriate initial orbital ele- 
ments of the lead satellite for several different time intervals. The 
software program developed to perform the numerical experiments is 
briefly described in Section 4.1. 

The constant parameters and the set of Initial orbital ele- 
ments of the trailing satellite which are held constant throughout the 
simulations are specified in Section 4.2. 
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The other sections of this chapter are devoted to the objec- 
tives of this investigation which are as follows: 

1. To establish a stable behavior in the relative motions under the 
influence of GEM-10B gravity field model (see Section 4.3). 

2. To determine the variations in the relative motion produced due to 
the effect of slight changes in the second degree zonal harmonic 
coefficient of the GEM-10B gravity field model (see Section 4.4). 

3. To determine the level of variations produced due to the effects of 
gravity field model error which is assumed to be the difference 
between the GEM-1 OB and the GEM-9 models (see Section 4.5.) 

The first objective is directly related to the investigation 
of the ’’stability” of the separation distance between the two satel- 
lites, wherein the GEM-1 OB gravity field model is assumed to represent 
the ’’best" approximate model representing the earth’s gravitational 
field. 


The second and the third objectives are investigated in con- 
nection with the fact that in the actual mission the satellites would 
sense the real effects of the earth's gravitational field which are 
expected to differ from the effects resulting from the a priori model. 
Therefore, it has been attempted to demonstrate the variations in the 
relative motions which would be ’’similar" to those in the actual mis- 


sion . 
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1 Software Description 

The software created to perform the numerical experiments of 
this investigation was designed to be implemented on the Control Data 
Corporation (CDC) Dual Cyber 170/750 computer system. The primary 
objective of the program was to perform the numerical optimization of 
the probl«n discussed in Chapter 3. Various software packages were also 
incorporated which were developed at the University of Texas at Austin. 
For example, the normalized harmonic coefficients of the geopotential 
expansion were processed by routines of UTOPIA (University of Texas 
Orbit Processor Incorporating Statistical Analysis) developed by Schutz 
and Tapley [1980], and the numerical integration of the equations of 
motion and the other differential equations required for the optimiza- 
tion algorithm were performed by using the DELIB software package writ- 
ten by McKenzie [19783. 

DELIB is a collection of six separate integration packages 
useful for the numerical solution of ordinary differential equations. 
Two of these integration packages are called ODE and ABFS, which are 
useful for solving a system of first-order differential equations. ODE 
uses the variable-mesh multistep method of Shampine and Gordon [19753, 
whereas ABFS uses fixed-order, fixed-mesh Adams method developed by 
Lundberg [19813. The results of the numerical integrations using the 
ODE package were more accurate compared with the use of the ABFS pack- 
age. However, the computation time efficiency was greatly improved by 
the use of the ABFS package due to the fixed-step Adams method with 


minimal overhead cost. In addition, both of the packages use the PECB 
(Predict Evaluate Correct Evaluate) technique where two function evalua- 
tions are made at each integration time step. Therefore, a remarkable 
computer execution time efficiency was also achieved by the use of 
psuedo function evaluations where only the two-body and the second 
degree zonal harmonic contributions to the accelerations of the satel- 
lites were reevaluated when the second function evaluation of the PECE 
technique was made. 

4.2 General Specifications 

In the numerical simulations, there were several parameters 
which, for simplicity, were held constant. Tliese parameters are; 

1. The epoch time (t^) and the angle of longitude of the Greenwich 

Meridian (a) at epoch which were set to be equal to zero. Hence, 

all of the simulations are generated with respect to, 

t s 0, and 
0 

a =0. 
o 

2. The constants which describe the physical and dynamical charac- 
teristics of the earth are specified in Table 4.1; wherein, is 

the equatorial radius, GM is the gravitational coefficient, and <a 

e 

is the mean magnitude of the earth’s rotation about the Z-axis. 


Initial Orbital Elements of the Trailing Satellite 


The particular set of initial orbital elements of the trailing 
satellite which are used for all of the numerical simulations are chosen 
according to the configuration vrtiich is shown in Figure 4.1. Tliis con- 
figuration suggests that the trailing satellite starts from the ascend- 
ing node of its orbit which is located at the altitude (h) of 160 kilom- 
eters on tlie X-axis. The initial circular speed assigned to the satel- 
lite is computed as follows 

r 

V = = 7808.034 m/sec. 

The respective initial orbital elements (Keplerian or the 
inertial components of the position and velocity) are computed which are 
given in Table 4.2, These initial orbital elements of the trailing 
satellite are held constant regardless of the various time intervals, 
for which the simulations are generated. 





Table 4.1 


Constants of the Earth 


GM 

(m^/sec^) = 

0.39860064 E-*-15 


(m) = 

6378145.0 

0) 

e 

(fad /sec) s 

0.72921151 E-04 


Table 4.2 

Initial Orbital Elements of the Trailing Satellite 
Used for all of the Numerical Simulations 


Inertial Rectangular 
Orbital Elements 

Keplerian 
Orbital Elements 

X (m) = 6538145.000 

a (m) = + 160000.0 

e 

y (m) = 0.0 

e = 0.0 

Z ( m) = 0.0 

i (deg) = 90.0 

X (m/s) = 0.0 

(deg) = 0.0 

u 

II 

o 

o 

(deg) r. 0.0 

Z (m/s) s 7808.034 

M (deg) r 0.0 


» Keplerian orbital elements are defined by; 

a = Semi-major axis 
e = Eccentricity 
i = Inclination 
w = Argument of perigee 
0, = Longitude of node 
M = Mean anomaly 
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4.3 General Behavior of the Relative Motions Characteristics 

The behavior of the relative motions between the two satel- 
lites is generally dependent on the initial orbital elements of the 
satellites. In this section, it is attempted to establish a stable 
behavior in the relative motions for a 12-day arc length, or time inter- 
val. All of the relative motion simulations are obtained with the use 
of the GEM-1 OB gravity field model which has harmonics complete to 
degree and order 36. The ntsnerical experiments which are performed to 
provide the simulations are as follows. 

A set of initial orbital elements for the lead satellite are 
chosen such that both of the satellites are located in the same orbital 
plane with equal circular speeds and radial distances from the center of 
the earth. Thus, having specified the initial orbital elements of the 
trailing satellite, the same initial orbital elements are used for the 
lead satellite except for a mean anomaly difference of 2.629 degrees 
ivhich gives approximately an initial separation distance of 300 kilome- 
ters. The respective inertial components of the position and velocity 
of the lead satellite as computed are given in Table 4.3. Using these 
initial set of orbital elements, the equations of motion for each of the 
satellites are numerically integrated to generate the simulations of the 
relative motions for a 1-day arc length Which are plotted as shown in 
Figures 4.2a through 4.2c (Case A. 1). Note that in Fig. 4.2a only the 
relative range variations from 300 km are shown. These plots clearly 
demonstrate the unstable behavior of the relative motions in the sense 
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of linearly increasing displacements with time. As can be seen in Fig- 
ures 4.2a through 4.2c there exists a 10 km decrease in the relative 
range, +1 m/s oscillations in the range-rate, and fluctuations in the 
relative acceleration magnitude which range from -0.002 to 0.001 m/s . 


Table 4.3 


Initial Orbital Elements of the Lead Satellite 
Used for the Case A. 1 


Inertial Rectangular 
Orbital Elements 

Ke Pierian 
Orbital Elements 

X (m) = 6531262.314 

a (m) = R„ + 160000.0 

e 

Y (ra) = 0.0 

e = 0.0 

Z (m) = 299921.040 

i (deg) = 90.0 

X (m/s) = -358. 174 

( 1 ) (deg) = 0.0 

Y (m/s) = 0.0 

Q (deg) = 0.0 

Z (m/s) = 7799.815 

M (deg) s 2.629 

1 
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The apparent periodic oscillations in the si-nulations can be 
explained as follows. The initially selected circular orbits change to 
elliptical orbits due to the perturbations of the gravity field, but 
assuming that the satellites follow each other in the same orbit, the 
lead satellite starts to slow down continuously after passing through 
perigee at a higher rate than the trailing satellite. Consequently, the 
separation distance decreases until the apogee location where it reaches 
the minimum value. However, an opposite behavior occurs after they pass 
the apogee location where the lead satellite starts speeding up faster 
than the trailing satellite, hence resulting an increase in the separa- 
tion distance until the next perigee location where the maximum separa- 
tion occurs. Thus, it can be concluded that during each 90 minute revo- 
lution of the satellites the separation distance goes through a cycle 
with a maximum and a minimum value corresponding to the perigee and apo- 
gee locations, respectively. 

Despite the large magnitude of the periodic oscillations, it 
is the secular change in the separation distance which essentially 
influences the stability of the system. The existing secular trend 
(Fig. 4.2a) is produced due bo the presence of time dependent nonlinear 
terms which would appear in a series expansion of the relative motion 
between the two satellites. Tlie effects of such terms can be eliminated 
by selecting initial orbital elements which tend to set the non-zero 
coefficients of the nonlinear terms to very small values or even ideally 
to zero. The selection of such initial orbital elements is accomplished 
by utilizing the optimization method which was developed in Chapter 3. 
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The procedure for this action is as follows. 

The initial guess for the orbital elements of the lead satel- 
lite required for the Newton-Raphson iteration method are chosen to be 
the same as those used for the previous simulations (Case A. 1). Then, 
following the steps of the estimation algorithm presented in Section 
3.3, several iterations are performed. The resulting initial state of 
the lead satellite at each iteration are summarized in Table 4.4, and 
the converged solution is given in Table 4.5. It should be noted that 
since the satellites were originally set in the inertial XZ orbital 
plane, the corrections made to the components of the position and velo- 
city of the lead satellite are mostly in the Y-direction. And, if the 
initial orbital elements had been chosen in the YZ plane, the correc- 
tions should have been made largely in the X-direction, which was vali- 
dated by other experiments. 

The significance of the resulting initial orbital elements is 
clearly reflected in the simulations which are generated again for a 1- 
day arc length as shown in Figures 4.3a through 4.3c (Case A. 2). These 
simulations, compared with the simulations of Case A.l, show that the 
secular trend in the relative range has completely been removed, and 
that a remarkable reduction in the amplitude of the oscillations has 
also been achieved, where now only +60 meters variations are present in 
the relative range simulation (Fig. 4.3a). 



55 


Table 4.4 


Results of the Iteration Process to Obtain the Appropriate 
Initial Orbital Elements of the Lead Satellite for 1-Day 
Arc Length with the Use of the GEM-1 OB Gravity Model 


Iter . 
No. 

X 

(m) 

■■ 

Z 

(m) 

lllllQl^^lllll 

MM 

|M| 

0 

6531262.314 

0.000 

299921.040 

-358. 174 

0.000 

7799.815 

1 

6531244.565 

-26134.483 

300101.156 

-358.839 

10.763 

7799.772 

B 

6531251.981 

-10383.628 

300209.453 

-358.962 

6.452 

7799.763 


6531253,126 

2051.539 

300214.241 

-358.966 

6.582 

7799.759 

■1 

6531251.485 

-631.542 

299977.100 

-358.687 

21.843 

7799.762 


6531252.199 

1154.880 

299967.175 

-358.674 

16.729 

7799.760 


6531252.237 

1778.375 

299960.775 

-358.667 

16.091 

7799.760 

B 

6531252.238 

1800.673 

299959.994 

-358.666 

16.100 

7799.760 

B 

6531252.238 

1800.544 

299960.014 

-358.666 

16.100 

7799.761 


Table 4.5 

Appropriate Initial Orbital Elements of the Lead Satellite 
Obtained for 1-Day Arc Length with the Use of the GEM-1 OB Gravity Model 
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In order to danonstrate the effects of the GEM-10B gravity 
field model for a longer time period, the obtained initial orbital ele- 
ments for Case A. 2 are used to generate the simulations for 3-day arc 
length which are shown in Figures 4.^a, b, c (Case A. 3). Obviously, the 
stable behavior is only present for the first day of the simulations and 
after that a divergent type behavior is developed where a 400 m decrease 
in the relative range is produced (Fig. 4.4c). However, using the same 
procedure as followed for Case A. 2, an appropriate initial state is 
obtained which provides stable behavior throughout the 3-day arc length. 
The initial state used for the iteration process of this case (Case A. 4) 
is chosen to be the one obtained for the 1-day aro length, and the 
resulting initial orbital elements are given in Table 4.6. The gen- 
erated simulations for this case are shown in Figures 4.5a, b, c. Note 
again the considerable reduction in the level of fluctuations which are 
achieved by the use of appropriate initial orbital elements. 


Table 4 .6 

Appropriate Initial Orbital Elements of the Lead Satellite 
Obtained for 3-Day Arc Length with the Use of the GEM-1 OB Gravity Model 


Inertial Rectangular 
Orbital Elements 

Keplerian 
Orbital Elements 

X (m) = 6531248.413 

a (m) = Rg + 159958.499 

Y (m) = 3635.284 

e = 0.540248 E-04 

Z (m) = 299944.245 

i (deg) = 89.884 

X Cm/s) = -358.630 

oj (deg) = 97.454 

Y (m/s) = 15.602 

n (deg) = 0.027 

Z (m/s) = 7799.768 

M (deg) = 265,181 
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The prediction of the relative motions characteristics was 
carried out for 6-day and 12-day arc lengths. Following identical pro- 
cedure as was done for Case A. 4, appropriate sets of ial orbital 

elements were obtained for 6-day and 12-day arc lengths which are given 
in Tables H.7 and 4.8, respectively. Tlie generated simulations for 6- 
day arc length using the obtained initial orbital elements for 3-day and 
6-day arc lengths are plotted as shown in Figures 4.6a, b, c (Case A. 5) 
and Figures 4.7a, b, c (Case A.6) respectively. Also, the generated 
simulations for the 12-day arc length using the obtained initial orbital 
elements for 6-day and 12-day arc lengths are shown in Figures 4.8a, b, 
c (Case A. 7) and Figures 4.9a, b, c (Case A. 8), respectively. Note the 
successive reduction in the levels of variations in the simulations for 
equal arc lengths. For example, comparing the simulations obtained for 
Case A. 8 with Case A. 7, a 50% reduction in the order of magnitude of the 
oscillations are achieved. 

The initial components of the position and velonity of the 
lead satellite obtained for various arc lengths are tabulated in Tables 
4.9a and 4.9b respectively. Note the consistent differences between 
each set of initial orbital elements. For longer t’me intervals the 
obtained initial states have a larger components for the position and a 
smaller components for the velocity in the inertial Y-direction. These 
initial orbital elements obtained for the lead satellite compared with 
the ones specified for the trailing satellite (Table 4.2) clearly indi- 
cate that the satellites move on very close orbits but not the same one. 






^.w r.- r-» - ...^ ^ ^ _ 




64 


Table 4 .7 

Appropriate Initial Orbital Elements of the Lead Satellite 
Obtained for 6-Day Arc Length with the Use of the GEM-10B Gravity ftodel 


Inertial Rectangular 
Orbital Elements 

Kepi erian 
Orbital Elements 

X (m) = 6531233.103 

a (m) = Rg + 159958.528 

Y (m) = 6495.894 

e = 0.485553 E-0 4 

Z (m) = 300078.439 

i (deg) s 89.901 

X (m/s) = -358.753 

oi (deg) s 96.741 

Y (m/s) s 13.198 

(deg) = 0.052 

Z (m/s) = 7799.775 

M (deg) = 265.895 


Table 4.8 

Appropriate Initial Orbital Elements of the Lead Satellite 
Obtained for 12-Day Arc Length with the Use of the GEM-10B Gravity Model 


Inertial Rectangular 
Orbital Elements 

Keplerian 
Orbital Elements 

X (m) = 6531193.906 

Y (ra) = 7170.572 

Z (ra) = 300266.230 

• 

X (m/s) = -358.912 

• 

Y (m/s) = 7.382 

• 

Z (m/s) s 7799.811 

a (m) = R + 159958.501 

© 

e = 0.401601 E-04 

i (deg) = 89.943 

0 ) ( deg) s 91 . 094 

n (deg) = 0.060 

M (deg) = 271.543 
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Figure M ,6a 

Relative Motions Characteristics for 6-Day Arc Length 
with the Initial Orbital Elements of Case A. 4 
Gravity Model: GEM-1 OB 
Case A. 5 
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Figure 4 ,Ba 

Relative Motions Characteristics for 12-Day Arc Length 
with the Initial Orbital Elements of Case A. 6 
Gravity Model; GEM-1 OB 
Case A. 7 

Relative Range Variations from 300 km (Ap ) vs. Time 
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Figure 4 .9a 


Relative Motions Characteristics for 12-Day Arc Length 
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Gravity Model: GEM-1 OB 
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Relative Range Variations from 300 km (Ap) vs. Time 
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Figure U .8c 
Case A, 7 (Continued) 
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Table 4 ,9a 


Appropriate Inertial Position Components of the Initial State 
of the Lead Satellite Obtained for Various Time Intervals 
with the use of the GEM-1 OB Gravity Model 


Arc Length 

X 

(m) 

Y 

Cm) 

Z 

(m) 

r 

(m) 

1-Day 
3-Day 
6 -Day 
12-Day 

6531252.238 

6531248.413 

6531233.103 

6531193.906 

1800.544 

3635.284 

6495.894 

7170.572 

299960.014 

299944.245 

300078.439 

300266.230 

+ 159991.971 
Rg + 159988.189 
R„ + 159981.270 
Re + 159951.441 


Table 4 .9b 

Appropriate Inertial Velocity Components of the Initial State 
of the Lead Satellite Obtained for Various Time Intervals 
with the use of the GEM-1 OB Gravity Model 


Arc Length 

' t 

X 

(m/s) 

• 

Y 

(m/s) 

• 

Z 

(m/s) 

V 

(m/s) 

1-Day 

-358.666 

16. 100 

7799.761 

7808.020 

3-Day 

-358.630 

15.602 

7799.768 

7808.024 

6-Day 

-358.753 

13.198 

7799.775 

7808.032 

12-Day 

-358.912 

7.382 

7799.811 

7808.068 
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of Error in the Second Degree Zonal Harinonic Coefficient 

The second degree zonal harmonic coefficient of the geopoten- 
tial expansion, given by Eq. (2.1), is denoted by J^, The motion of a 
satellite is dominantly perturbed by the effect of the Jg which 
represents the earth’s oblateness. Therefore, in this section the Jg 
coefficient of the GEM-10B gravity model is purposely altered by a 
slight amount in order to demonstrate the produced effects in the rela- 
tive motion simulations due to this change. The simulations procedures 
are as follows . 

The normalized coefficient of the GEM-1 OB gravity model has 
the value of -0.4841655 E-03 and was increased by 1.0 E-08 for the simu- 
lation to produce the effect of an error in the coefficient. The ini- 
tial orbital elements of the lead satellite, which were obtained for 
Case A. 2, were used to simulate the relative motions. The difference 
between the obtained simulations for this case (Case B. 1) and the ones 
obtained for Case A. 2 are shown in Figures 4.10a through 4.10c. 

The same initial orbital elements were used to obtain an 
appropriate set of initial orbital elements for the lead satellite which 
remove the small secular trend in the separation distance which is 
apparent in Fig, 4.10a. This was accomplished by applying the same 
optimization method as before. The obtained initial state of the lead 
satellite for this case (Case B.2) is given in Table 4.10. Once again, 
the simulations obtained for this case are differenced with the simula- 
tions of Case A. 2 and the results are plotted as shown in Figures 4.11a 


through 4,11c. These simulations compared with the simulations of Case 
B.1 indicate a well-behaved periodic oscillations which represent the 
significance of the obtained initial orbital elements. 


Table 4.10 


Appropriate Initial Orbital Elements of the Lead Satellite 
Obtained for 1-Day Arc Length with the use of 
Gravity Model; GEM-1 OB +AJg 


Inertial Rectangular 
Orbital Elements 

Keplerian 
Orbital Elements 

X (m) = 6531252.184 

a (m) = R^ + 159959. 107 

Y (m) = 1813.933 

e = 0.568937 E-04 

Z (m) = 299960.613 

i (deg) = 89.881 

X (m/s) = -358.667 

0) (deg) = 97.698 

Y (m/s) = 16.080 

(deg) = 0.010 

Z (m/s) s 7799.761 

M (deg) = 264.938 
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the Initial Orbital Elements of Case A. 3 
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F igure 4.11a 

Relative Motions Perturbations Due to aj *10"^ 
of the GEM-10B Model for 1-Day Arc Ler^th 
with Appropriate Initial Orbital Elements 
Case B.2 
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4.5 Effects of the Gravity Model Error 

So far in this chapter it has been assumed that the GEM-1 OB 
gravity model represents the most accurate model of the earth’s gravita- 
tional field and all of the relative motions simulations were obtained 
with the use of GEM-10B model. In addition, a special case was also 
treated in Section 3.4 where a set of simulations were provided in order 
to demonstrate the effects of a possible error in the coefficient of 

the GEM-10B model (Cases B.1-2). However, in this section the effects 
of a full gravity model error in the relative motions are demonstrated 

in which the GEM-9 model which has harmonics complete to degree and 
order 22 is assumed to represent the "true” gravity model of the earth. 

The initial orbital elements of the lead satellite which were 
obtained for Case A. 2 (given in Table 4.5) are used to simulate the 
relative motions under the influence of the GEM-9 gravity model. The 
resulting simulations are plotted as shown in Figures 4.12a through 
4.12c (Case C.1) and the differences between these simulations and the 
ones obtained for Case A. 2 are shown in Figures 4.13a through 4,13c 
(Case C.2). Note that because of the different energy of the GEM-9 
model, the stable behavior of the relative motions in Case A. 2 is com- 
pletely changed to an irregularly secular behavior in which a 600 meters 
error in the separation distance is produced (Fig. 4.13a). 
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Figure 4 ,12a 

Relative ^tation3 Characteristics for 1-Day Arc Length 
with the Initial Orbital Elements of Case A. 2 
Gravity Model; GEM-9 
Case C. 1 



Figure 4.13a 


Relative Motions Differences between GEM-10B and GEM-9 Models 
for 1-Day Arc Length with the Initial Orbital 
Elements of Case A. 2 Used for each Model 
Case C.2 
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F Igure 4 ,12c 
Case C.1 (Continued) 


Relative Acceleration Magnitude (p) vs. Time 
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Figure 4 ,13c 
Case C.2 (Continued) 

• • 

Relative Acceleration Magnitude Differences (5 0) vs. Time 



TIME PAST EPOCH (DAY) 




80 


The same initial orbital elements are used to obtain appropri- 
ate set of initial orbital elements such that the generated simulations 
best fit the simulations of Case A.2 in a least squares sense. This is 
accomplished by minimizing the performance index which can be expressed 
as 


1-day 2 

^ ^ ^ ^GEM-IOB “ ^GEM-9^ 


where represents the relative range between the two satellites 

computed at some instant of time with the use of GEM-9 model; whereas, 
*^GEM 10B corresponding relative range computed at the 

same instant of time with the use of GEM-10B gravity model . The pro- 
cedure to minimize the above function is basically the same as that 
presented in Section 3.3. Therefore, after performing several itera- 
tions, the final adjusted initial orbital elements for the lead satel- 
lite are obtained vdiich are given in Table 4.11. Then, the obtained 
initijii orbital elements are used to generate the simulations with the 
use of GEM-9 model. These simulations as well as their differences with 
the simulations obtained with the use of the GEM-1 OB model are shown in 
Figures 4.14a through 4.14c (Case C. 3)* and Figures 4.15a through 4.15c 
(Case C.4) respectively. Note the obtained +16 meters error in Fig, 
4.15a as compared with the previously existed 600 meters error in Fig. 
4.13a. 
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Table 4.11 

Appropriate Initial Orbital Elements of the Lead Satellite 
Obtained for 1-Day Arc Length with the use of 
Gravity Model: GEM-9 


Inertial Rectangular 
Orbital Elements 

Keplerian 
Orbital Elements 

X (m) = 6531251.338 

a (m) = Rg + 159962.234 

y (m) s 1826.607 

e = 0.57^876 E-04 

Z (m) = 299945.957 

i (deg) = 89.881 

X (m/s) = -358.655 

0 ) (deg) = 96.936 

Y (m/s) 2 16.081 

n (deg) = 0.011 

Z (m/s) = 7799,765 

M (deg) = 265.700 


Having obtained the appropriate initial orbital elements, the 
error in the relative motions is predicted for six days. The simula- 
tions for a 6-day arc length using the GEM-9 model are shown in Figures 
4.16a through 4.16c (Case C.5), and their differences with the ones 
obtained with the use of the GEM-1 OB model are shown in Figures 4.17a 
through 4.17c (Case C.6). Note again the 600 meters error in the 
separation distance which is produced during the additional five days of 


the simulations. 
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Figure U ,14a 

Relative Motions Characteristios for 1-Day Arc Length 
with Appropriate Initial Orbital Elements 
Gravity Model: GEM-9 
Case C.3 


Relative Range Variations from 300 km (Ap ) vs. Time 
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Relative Motions Differences between GEM-10B and GEM-9 Models 
for 1-Day Arc Length with Appropriate Initial 
Orbital Elements Used for each Model 
Case C.4 
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Case C.3 (Continued) 
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Figure ^.16a 

Prediction of the Relative Motions Characteristics for 6-Day 
Arc Length with the Initial Orbital Elements of Case C.3 
Gravity Model: GEM-9 
Case C.5 


Relative Range Variations from 300 kra (Ap) vs. Time 



Figure ^ ,17a 

Prediction of the Relative Motions Differences between 
GEM-10B and GEM-9 Models for 6-Day Arc Length with 
the Initial Orbital Elements of Case G.4 
Case C.6 
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CHAPTER 5 


CONCLUSIONS AND RECOMMENDATIONS 

5. 1 Summary 

This preliminary study was oonduoted by performing mmerous 
numerical experiments which provided the simulations of the relative 
motion characteristics between the two low altitude satellites. The 
simulations demonstrated the effects of two complete GEM-10B and GEM-9 
gravity field models. Furthermore, there were various intervals of time 
(arc lengths) for which the simulations were generated. 

Throughout the numerical experiments, the method which was 
developed in Chapter 3 was used to determine appropriate set??, ^?f initial 
orbital elements for the lead satellite whereas the initial orbital ele- 
ments of the trailing satellite were held constant. The initial orbital 
elements of the lead satellite were obtained such that the two satel- 
lites would sense the same during the itntlre length of the simu- 
lation time. Consequently, in the generated simulations, only the 
rainimun level of variations of the relative motions were demonstrated. 
In thi,? manner, the ability of the method was well acknowledged in 
Chapter 4 wherein a maximum of 12-day simulation time was considered . 
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5 . 2 Conclusions 

From the analysis of the numerical simulations provided in 
this investigation, the following conclusions can be made. 

The simulations provided in Tection 4,3 (Cases A. 1 through 
A, 8) demonstrated the fact that the stability of the system of two 
satellites is very sensitive to the initial orbital elements of the 
satellites. In the various cases which were considered, the stable 
behaviors were only possible to obtain with the appropriate selection of 
certain accurate set of initial orbital elements corresponding to the 
desired intervals of time. Moreover, it was also shown that the initial 
orbital elements obtained for an specified interval of time would ncl 
necessarily promise an stable behavior for a longer simulation time. In 
other words, an appropriate set of initial orbital elements were 
required in order to establish stability for the additional portion of 
the simulation time. Therefore, it can be concluded that in the actual 
mission it may be very difficult to hold the separation distance to 
small fluctuations; consequently, periodic adjustments to the orbits of 
the satellites may be required . 

In addition, the simulations provided in Section 4.4 and 4.5 
demonstrated the level of perturbations in the relative motions which 
were produced due to the effects of gravity field model error. Of 

course the level perturbations were also greatly reduced by the use of 
appropriate set of initial orbital elements which counteracted the 
effects of the error in the gravity field model. 
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Given that the GRAVSAT mission is expected to provide an accu- 
rate recovery of the geopotential model and that in the actual mission 
precise initial orbits for the satellites can not be practically 
achieved, it is essential that in the derivation of accurate observation 
equations (required for processing the actual data) which determine the 
quality and precision of the results, a considerable caution should be 
taken whenever a simplifying assumption is used or an approximation is 
made. 

5.3 Recommendations 

Tlie following recommendations are suggested in order tn 
improve the precision of the computed relative motions characteristics 
and the realism of the simulations. 

The level of variations of the relative motions can be com- 
puted more accurately by converting the current single precision coded 
software program into double precision. 

The only forces acting on the satellites considered in this 
study have been the central body gravitational forces. However, it is 
of importance to include the perturbations which will be produced due to 
other forces such as the effects of polar motion, direct and reflected 
solar radiation pressure and other non-gravitational forces. However, 
the drag compensation mechanism proposed for GRAVSAT should accommodate 
all surface forces, including radiation pressure. It should be noted 
that by including any of the above improvements the computer execution 
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time of the already considerable time consuming process will be 
increased. 


It is suspected that the simulations of the components of the 
relative motions in the radial, transverse or along-track, and normal or 
cross-track directions can be more descriptive of the relative motion 
characteristics as compared with the simulations along the line of sight 
between the two satellites. 

Finally, as mentioned previously, the simulations provided in 
this study were obtained with respect to an arbitrary set of initial 
orbital elements for the trailing satellite which, regardless of the 
various intervals of time, were held constant. However, it is suggested 
that for further studies these initial orbital elements be selectea such 
that for any given interval of time and gravity field model, when the 
orbit of the trailing satellite is determined, it best fits the one in 
which only the two-body forces are modeled. This can be accomplished by 
using the same minimization method except that the constant separation 
distance should be set to zero. It is expected that the simulations of 
the relative motions obtained with the use of such initial orbital ele- 
ments will demonstrate anoother behavior in the relative motions and 
that the fluctuations will be distributed more evenly along the simula- 


tion time. 
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